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Abstract 


This report describes the development and applications of multiblock/ 
multizone and adaptive grid methodologies for solving the three-dimensional 
simplified Navier-Stokes equations. The program was initiated in 1987 focusing on 
developing a three-dimensional plume code to simulate the aerodynamic 
characteristics of a jet. issuing from nonaxisymmetric nozzles. Previously, Abdol- 
Hamid et. al. introduced the single zone version of the present code (PAB3D-vl) where 
the parabolized and simplified Navier-Stokes equations were solved. The code was 
tested and compared with the experimental data for axisymmetric underexpanded and 
overexpanded supersonic jet flows and transonic flow around a nonaxisymmetric 
afterbody. 

In the present report, adaptive grid and multiblock/multizone approaches are 
introduced and applied to external and internal flow problems. These new 
implementations Increase the capabilities and flexibility of the PAB3D code in solving 
flow problems associated with complex geometry. 



1. Introduction 


A single block solver can be used efficiently to simulate simple aerodynamic 
configurations. Among various methods offered by many researchers, Abdol-Hamid 
1,2,3 introduced the single block version of PAB3D code to simulate underexpanded and 
overexpanded supersonic jets issued from round and rectangular nozzles. Abdol-Hamid 
and Compton^ used the PAB3D code to simulate external flow around a nonaxisymmetric 
nozzle at a Mach number of 1.2. Pao and Abdol-Hamid^ used the single block with 
adaptive grid to simulate underexpanded supersonic jet flows issued from round, 
square, and elliptic nozzles. 

As better computational methods and powerful computers are available in 
recent years, computational fluid dynamics (CFD) has become one of the important 
tools in improving aircraft design (6,7). Until recently, the use of CFD was limited to 
simple geometries. Future aircraft (fighter or transport) will have very complex 
geometries and are difficult to handle with a single zone structured grid. Either 
unstructured or multiblock/multizone structured grids are attractive approaches for 
solving viscous flow problems with complex configurations. Even though the 
unstructured grid is much easier to generate, it requires more computational time and 
memory for solving the Navier-Stokes equations per grid point. With the capability of 
the supercomputers of today, the multiblock/multizone approach is a flexible method 
which can handle very complex configurations. 

The advantages in using the multiblock/multizone approach are: 

1. Simple grid generation for complex configurations. 

2. Flexibility to use a different CFD approach for each block: 

a. Numerical technique (space marching algorithms for supersonic flows 
and time-dependent algorithms for subsonic and separated flows). 

b. Different topology for each block (polar, cartesian, etc.). 
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c. Adaptive grid in regions where the dependent variables and their 

gradients change their strength and location. 

3. Less memory as each zone is solved independently with appropriate 

boundary conditions. 

This report describes the capabilities of an improved version of the PAB 3D-vl code 
reported in references 1 to 4. This improved code, named PAB 3D-v2, includes options 
for three different numerical schemes to solve the simplified Navier-Stokes equations. 
The three schemes are: the flux-vector-splitting scheme of van Leer®, the flux- 
difference-splitting scheme of Roe® and a modified Roe scheme (space marching 
scheme)^*®. Four different turbulence model options are also included in PAB 3D-v2. 
The first of the four, the Baldwin-Lomax 1 ® model, is a two-layer algebraic model which 
follows the pattern adopted by Cebeci 1 1 but avoids the necessity of determining the 
boundary layer thickness. The second, the Johnson and King model ^ a s extended to 
three-dimensional flows by Abid*® and Abid el. al.^. is a two-layer hybrid eddy- 
viscosity Reynolds shear-stress model in which a simplified ordinary differential 
equation for the maximum Reynolds shear-stress is solved. The third, the Goldberg 
model as modified by Goldberg and Chakravarthy can be considered as a three- 
layer turbulence model where the third layer is used to simulate the separated regions of 
the flow. The last is the mixing length turbulence model^ with the option of including a 
compressibility correction factor introduced by Cheuch^ ^ Two different external and 
one internal flow problems are used to test the various code capabilities. 

One important problem for CFD applications is the prediction of the shock-cell 
structure of underexpanded and overexpanded supersonic jet flows. Understanding the 
effect of shock-cell structure and interaction of a supersonic jet with the external 
stream is essential for the design of future aircraft. Also, the nozzle exit geometry plays 
an important role in designing fighter aircraft for maximum maneuverability over a 
wide range of Mach numbers 1 ®*^ Developing an efficient computational technique is 
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Important to fully understand the flow characteristics of these nozzles. At the present 
time, there are few codes available to predict the aerodynamics of three-dimensional 
shock containing jets. Wolf et. al. developed a three-dimensional code (SCIP3D 2 3) for 
analyzing the propulsive jet mixing problem. Anderson and Barber 24 also developed a 
three-dimensional Parabolized Navier-Stokes procedure for calculating the heated 
subsonic and supersonic jet. This code was used to simulate the jet mixing rate for 
axisymmetric, rectangular and splayed nozzles operated at design conditions. Abdol- 
Hamid 2, 3- 4 introduced a space marching scheme, which is based on modifying the Roe's 
scheme, to get an accurate solution to the simplified Navier-Stokes equations for 
supersonic flows with a single time sweep. This scheme was successfully used to 
simulate underexpanded supersonic round and square jet flow problems 2 ’^. Pao and 
Abdol-Hamid^ introduced a new adaptive grid for analyzing the aerodynamic of shock- 
containing single jets. They used this technique to simulate round, square, and elliptic 
jet flows. The adaptive grid is used to accurately describe the shear layer and detect and 
track the movement of the shock system for underexpanded supersonic jets. In the 
present report, adaptive grid and multiblock capabilities included in PAB 3D-v2 are 
utilized to simulate round, square, and elliptic supersonic jet flows. 

Another group of underexpanded supersonic jet flow which involving the 
Internal and external flow regions for a special family of jet nozzle is analyzed In this 
report. These examples are designed for showing the flexibility of the PAB3D-v2 code in 
handling mixed boundary conditions over a block interface. The nozzle configuration 
can be described as a circular pipe section followed by five equally spaced tabs. Each tab 
is simply the extension of an arc segment of the circular pipe for a certain length in the 
downstream direction. Each arc segment, representing the width of the tab, is 1/10 of 
the full circle. For this family of configurations, only two grid blocks are needed for 
calculations using the PAB3D-v2 code. It is estimated that at least 30 percent of 
computer resources are saved by such structural simplicity when compared to typical 
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multiblock codes. Results of analysis using PAB3D-v2 for these nozzles are 
qualitatively similar to the experimental results obtained by Wlezien et al 44 for 
nozzles with 1. 2, 4 and 8 tabs. In general, the results show that the tab nozzle 
configuration allows rapidly establishment of a pressure equilibrium between the 
underexpanded jet flow and the ambient free stream. The jet plume is found to have a 
higher spreading rate and a lower core flow Mach number as compared to a similarly 
underexpanded supersonic jet issuing from a circular nozzle without tabs. 

Finally, PAB3D-v2 was used to predict (he aerodynamics of an afterbody at 
transonic speed. In fighter development programs, a great amount of effort is spent in 
analyzing the afterbody flowfield to efficiently integrate the nozzle and airframe. For 
analyzing this complex flowfield, computational fluid dynamics is becoming 
increasingly useful. Previous applications of computational fluid dynamics to the 
afterbody problem include numerical techniques ranging from panel methods to 
Navier-Stokes solvers^5-33 Abdol-Hamid and Compton 4 used four different 
numerical algorithms and three different turbulence models to solve the three- 
dimensional Navier-Stokes equations for supersonic flow over a nonaxisymmetric 
nozzle. Three of the algorithms were contained in the PAB3D-vl arid PAB3D-v2^'^ and 
the other in the CFL3D code^^>^ 4 '^®. In the present report, the mulliblock/multizone 
approach in PAB 3D-v2 is utilized to simulate the flow over this nonaxisymmetric 
nozzle at a Mach number of 0.8 using a coarse grid. Also, the performance of the three 
turbulence models using a fine grid topology in simulating supersonic flow are 
compared with experimental data. 

2. Governing Equations 

The governing equations under consideration here are the Reynolds-averaged 
Simplified Navier-Stokes equations obtained by neglecting all streamwise derivatives, 
3/3J;, of the viscous terms. The resulting simplified Navier-Stokes equations are 
written in generalized coordinates and conservation form as 
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where. 



E = j($ t Q + $,E + $ y F + $ 8 G) 


F -j(^tQ^ x E + H y F + n z G) . I(ti x E v +Tl y F v + Ti z G v ) 
^ = j(^t9 + ^x E + Cy F +C Z G) - — (C X E V +CyF v +C Z G V ) 


( 2 ) 
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pu 

pu 
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pv 

, E= puv 

pw 

puw 

e 

(e + P)u 


pv 

pw 

puv 

puw 

pv2 + P 

, G = pvw 

pvw 

pw^ + p 

(e + P)v 

(e + P)w 


( 3 ) 


In these equations, p is the density, u, v, and w are the components of the velocity 
in the x, y, z directions, respectively, and e is the total internal energy per unit volume. 
The pressure, P, is related to the energy by 

P =(Y-l)[e-|-p(u 2 +v 2 + w 2 )] (4} 

o 

*xx 
Ev = "txy 
*xz 

Y rYT 

— K ~ + UT x „ + VT xv + WX 
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Gv = tyy (5) 
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where 


^ = ^(x,y,z,t) = Streamwlse (marching) direction 
T| = Ti(x,y,z,t) = Normal direction 
£ = £(x,y,z,t) = Spanwise or circumferential direction 
J is the Jacobian of the transformation given by 


j =a(fn.Q/c(x.y,z) 


_ _ .. I o 3u dv dw 
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dx dy dz 
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_ * I 2 _dv _ 3u _ dw 
° ^ dy dx ~dz 

dw du dv 
dz dx dy 


— — ill O u w _ wu v 

ZZ “ g ^ I Z ^ ^ ~ -V 


p = p L +W 
M- = WL 



dy 

dx 



afterbody calculations 
jet and nozzle calculations 


(6) 


PL and pt are the laminar and turbulent viscosity respectively. In the present 
investigation, the turbulent viscosity is evaluated using two algebraic turbulence 
models which are described subsequently. 
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The Parabolized Navier-Stokes (PNS) equations are obtained from the 
governing equations when the unsteady terms are omitted and the following 
assumptions are enforced: 

1. The streamwise velocity component is everywhere greater 
than zero. 

2. The pressure gradient term in the streamwise direction 3P/ 
is either omitted or treated with other techniques to avoid a 
complex eigenvalue. 

In the present investigation, the technique of Vigneron et. al.^ 7 is adopted to 
suppress the departure solutions associated with the elliptic behavior of the equations. 
Vigneron et. al.^ 7 show that PNS equations are hyperbolic-parabolic provided that the 
streamwise convected flux vector is replaced by 

E=[pU,puU + S x coP.pvU + S y ioP,pwU + S z cop,(e +P)u] ^ 

where 

(0=1 M- > 1 

2 

oyMt 

= L - L ~r M 4 <1 

i+(y-i)m* 

and, a is a safety factor to account for the nonlinearity of the governing equations. A 
value of 0.95 is used in the present calculations. 

3. Turbulence Models 

The Baldwin-Lomax, Johnson-King, and Goldberg turbulence models (for wall 
boundary problems) and mixing length turbulence models (for shear (low problems) are 
briefly described in this section. 
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3.1 Baldwln-Lomax Model 


The Baldwin-Lomax^O model Is an algebraic two-layer turbulence model which 

follows the pattern adopted by Cebeci 1 1 . The model Is developed for thin-layer, two- 

and three-dimensional calculations. The turbulent viscosity is evaluated as follows: 

Mt = Mi q * Tlcros 

Mt = (Mt)o ti > Tlcros 

where q is the normal distance from the afterbody surface and q cros is the smallest 
value of T| at which values of and (pi) 0 are equal. For the inner-layer: 

(pt)i = pt 2 M (8) 

where, i = kq [l -exp (- q + /A + )] 
and, k = 0.4, A + = 26 

+ _ Pw^W 1 ! 

^ Hw (9) 

where x w = wall shear stress 

l M l is the magnitude of the vorticity. 


There are two different ways to calculate I 

/- V . - X ... Y 4 . \2 


M3D 


3u 3v f 3v dw ( dw du 

dy dx ) + ^dz dy J + (,dx dz 


and the thin-layer approximation 


Nt-l = 


du 

T1y a^ _T1x 


dv Y [ dv dw 

dq J + [^ z dq T * y dq 


to i ; the three-dimensional form: 

1/2 

(10a) 


f dw du 


\ 2 


1/2 


(10b) 


For the outer-layer: 

(M^o = KC cp p f'wake ^'kleb M 
where, C^p = 1.6, k = 0.018 

F wa ke = the smaller of 


^max ^max 
^-wk^maxU / ^max 


( 11 ) 


U = u^ +y2 + w^» C w k= 0.25 

q max is the q location corresponded to the maximum value (F max ) where F is 
calculated by 
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( 12 ) 


F ( Tj) = T||co| [l - exp(- ti + /A + )] 
and 



(c > 

6 ** 

1+5.5 

klcb 


^Imax j 



_ 


_ 


Ckleb - 0.3 

3.2 Johnson-King Model 


The Johnson-King model 12 ’ ^ j s a two-layer hybrid eddy-viscosity/Reynolds- 
shear-stress model. A simplified ordinary differential equation for the maximum 
Reynolds shear stress (x m ) is used to determine the change in the turbulent viscosity in 

the streamwise direction. The initial values of the Reynolds shear stress (along each 
line normal to the afterbody) have to be evaluated with some other approaches. The 
original model was developed for 2-D flows in which a 1-D equation for x m is to be 

solved. Abid 1 ® and Abld et. al. ^ extended the Johnson-King model 12 for the 
application to 3-D flows. In the Johnson-King model, the turbulent viscosity is 
expressed as 


1/2 

(Pt)i=fcn*m 


+ + 
l-exp| -q, /A, 


(15) 


where 
+ 1/2 
*ll = P^l^max / Pi 

Aj + =17 


T max = maximum Reynolds shear stress/density 
t = pjwl/p 


(16) 
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The outer eddy viscosity Is the same as the one used for the Baldwin-Lomax model 
(equation 1 1) but multiplied by a correction factor a. However, k takes a value of 0.0168 
as suggested by Abid et. al. 14 . The a factor provides a link between the eddy viscosity 
evaluated by equation (16) and x m . x m is evaluated by solving the 2-D ordinary 
differential equation, which can be written in the following finite volume form: 


W 

¥ + ufV + r=0 

where 


(17) 


"*‘ # f 2 W m 

3 i f g "j C P L ml 1 ~P 1/2 | Vq I 

2 L mV^eq J 8aj (0. 7 - q m /5) U m 


(18) 


where aj = 0.25, Cp = 0.5 
Lm = min (0.4 q m . 0.098) 

First, the Baldwin-Lomax model is used to supply the initial values for tni at each 
streamwise location, and a is set to 1. Then, at the following time steps, equation (17) is 
solved for t m using an upwind-scheme, and o is updated as follows 


a t + At = o t x in p niax /(p M) 

v ’ max ( 19 ) 

In the region where a is less than unity, the value of |l - o 1/2 l (equation 14) is set to zero. 
3.3. Goldberg Modification 

Goldberg *5 and Goldberg and Chakravarthy 1 ® introduced a modification for 
boundary layer turbulence models, which is designed to simulate the separation bubble 
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in the flow. The modification consists of a simple formula for the distribution of the 
eddy viscosity within the separation bubble. The Baldwin-Lomax model is used outside 
the separation bubble and the edge of the separation region "t| s " is treated as a wall 
boundary. Thus, the q in equations 9, 11, 12 and 13 is replaced by ” 11 -%". 

For high Reynolds number flows, Goldberg and Chakravarthy^® show that the 
viscosity in the separation regions can be evaluated as 


(Pt ) s = C]U S % V^[A(q / % )+ B]VG 

and 

(fit) v = c 2 u s T lsVPwP q s < q < q cr 

where 

G = [1 - e^fiAW 2 ] / [1 - e A 


o<ti<ti s 


A = — 


, 9/5 

( \ 

C^/2j , B = 

2 


J 


-A 


( 20 ) 


( 21 ) 


,, _ l U t.m^max 

8 v p 

Pt.m = Pt l a> max 

C x = 0.353, C 2 = o. 188, + = 0.5, C (l = 0. 7 

In the separation layer, equation (16) is used up to Tier . the smallest value of q at which 
values from equation (20) and Baldwin-Lomax model are equal. Beyond q C r. (he 
Baldwin-Lomax model is used. 

3.4 Mixing Length Model (ML) 

This is an algebraic eddy-viscosity turbulence model which is based on the 
Prandtl hypothesis. The turbulence viscosity is evaluated as 



where i is the turbulence length scale, 
i = 0.11tig 
nG = T|2-m 
where at t|i 


U-U«o 

Uo-U„ 


= 0 . 1 


and at T |2 


(23a) 


U-U, 

U 0 -U 


= 0.9 


(23b) 


where U = £xu + £yv + ^w 

is the contravarient velocity component in the streamwise directions, U „ is the 
external flow velocity and U 0 is the jet centerline velocity. 

3.5 Compressibility-Corrected Mixing Length Model (ML-CC) . 

It is well known that turbulent mixing rales are reduced for supersonic flows in 
comparison to subsonic flows. Chuech et. al. 17 introduced a compressibility 
corrections factor and used it to modify the turbulent viscosity: 
p t = KpiM (24) 


where k 

= 1.0 M c < 0.55 

= 2.03 - 1.87M C for 0.55 ^ M c ^ 0.95 

= 0.25 Me > 0.95 

where M c is the connective Mach number of the mixing layer. 


4. Computational Methods 

The three computational schemes presented in this report are basically 
implicit, upwind, and constructed using a finite volume method. The diffusion terms 


are centrally differenced and the inviscid flux terms are upwind differenced in these 
schemes. Associating the subscripts i,k,j with £, r|, £ directions, a numerical 
approximation to Eq. ( 1) may be written in the following form: 


/ « \ A uti iit i a jin ^ jit i iiti ^ n + 1 

(Qi.k.j), + E j — E i +F j — F ] +G i — G j =0 

' h i+ 2- k ’J * ,k+ 2’1 1k -2’ J ‘• k -J + £- lk 'j-5- 

The fluxes at (n + 1 time iteration) are linearized as 
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^ n + 1 n n 

G =G + 3Q aQ 


Then, equation (25) is written as. 
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In the present code, two flux-splitting schemes are used to construct the convective flux 
terms in equation (26). 

The variables Q + , Q" are defined by an upwind biased one parameter 

family 


9 , + 1 k . " 9 <*.j + 4 {I 1 ~ + 1 1 + 9 i.k.j 

2’ J 

(27) 

9 , + I. k>J = 9 J + lk -J “ 4 {( 1 + k ^) + ( 1 _k ^) A 4 9l - k -J 

These variables can be either the conservative or primitive variables. Also, Q+ and Q- 
represent the right and left variables, respectively, in reference to the cell face, 
where 

A$Qi.k,J = Ql + l.k.j — Q l.k.j ’ A !-Q l.k.j =9 l.k.j- Qi - l.k.J 

<)> = 0 first order fully upwind 


k^ = -1 
<t>=l 


second order fully upwind 



♦ =1 


third order biased upwind 


However, to ensure monotonic interpolation for the third order interpolation in the 
vicinity of a shock, a min mod limiter is used as follows: 

VQ = min mod (VQ.bAQ) 

(28) 


AQ = min mod (AQ.bVQ) 


3 - k £ 

where b is a compression parameter, b = - 
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It should be mentioned that the splitting procedures are only used for the 
inviscid convection parts of the flux vectors (f and 6) . A second order, central 

difference is used to represent the diffusion (viscous) terms. 
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Jj , J 2 , J 3 , Vj and V 2 are completely described in Ref (1). 

The implicit upwind /relaxation algorithm of Newsome et. al.®® is used to solve 
the governing equations. This can be achieved through a series of alternative sweeps in 

n n + 1 

the streamwise direction. For a forward sweep. AQ, _ , k j is known and AQ 4 + lkj is set 

n n+1 

to zero. For a backward sweep. AQ, + , k j is known and AQj _ lkJ is set to zero. 

Finally, equations (29) are approximately factored and can be written in the 
following compact form: 


L + 8 r 


f 9F | 3F V ^ ~' { 

dQ + 3Q j 


L + 8, 
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n + 1 

AQ = R.H.S 


(30) 


where, 

L=™ + B, 
At 5 


4.1 van Leer flux vector-splitting scheme 


The first of the two schemes is the van Leer® flux-vector-splitting method. The inviscid 
terms of the flux vectors (e, F and 6) are split according to their contravariant Mach 

number. 

U V W 

(M^.M^andM^), defined as = - and M ? =— ( 

where 

U = (S x u + S y v + S z w) / S 
V = (R x u + R y v + R z w) / R 

W = (T x u + T y v + T z w) / T 

2 2 2 2 
S = S X + Sy + S z 

2222 2222 
R = R x + Ry + R z and T = T x + T y + T z 
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As an example, for supersonic flows in the x direction 
® and E = 0. for 


^ _ 0. and E ~ E. for 
and for subsonic flows, -1 < < 1 


± 

E = S 


EL,{s x (-U±2a)/YS + u} 
E* asa {s x (-U±2a)l7S + v} 

Emass {s z (~u ± 2a)/ yS + w} 


E, 


energy 
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< - 1 


( 32 ) 


where 


E m ass = ±Pa ( M 5 ± l) 2 / 4 
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4.2 Roe's flux-difference -splitting scheme 

The second scheme is the Roe's flux-difference-splitling method^, which solves the 
approximate Riemann problem. For example, the interface flux in the streamwise 


direction is evaluated as. 



E +i = I[E(Q L ) + E(9 R HAI(9 R -9 L )] t + i 
2 2 


( 34 ) 


where 9 l (9 ) and 9 r (Q + ) are either primitive or conservative variables to the left and 

the right of the cell faces, and A is the Roe-averaged flux Jacobian matrix: 

The last term in equation (34) I Al (9 r-9l) Is defined as: 


a 4 


iai(9 r -9l) = 


ua 4 + k x a 5 + a 6 
ua 4 + k y a 5 + a 7 


wa 4 + k z a 5 + a 8 

.2 

Ha 4 + a 5 u + a 8 u + a 7 v + a 8 w - 


(35) 


u = 


u L*JpL + u R ^ 




(36) 


Also, p.v.w and H are evaluated using formulas similar to eq. (36) and 


c = (Y - 1)H - 


2 2 2 
U + V + W 
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where 
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a 7 = Kl ul p (Av -k y Au) 
a fl = Klulp(Aw - k z Au) 
and 


K = 


12 2 2 

yk x + ky + k z 
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For fully supersonic flow, the Information (disturbance) can only travel in the flow 
direction according to equation (34) 


^i+l/2,k.J 



Jl 


+ 1/2. k.J 
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(39) 


2 ( Ay A 

°^ vc "' = Backward difference 

where is the average Mach number evaluated using the Roe-averaged method. The 

same result is achieved by setting the state variables at the i + 1 and i + 2 planes to be 
equal to the ones at the J-plane, 

9i + 2 = Qi+l=Qi fy>l (40) 

For a large number of supersonic/subsonic mixing flow problems, the 
downstream influence (contributed from the subsonic regions) can be neglected when 
compared with the upstream effect. For these cases, equation (18) can even be used in 
the subsonic regions to obtain an accurate solution. 

With this approach, a simple modification to the Navier-Stokes solvers can be 
made and a wide range of problems can be simulated with a single sweep. This approach 
does not require any changes in the governing equations restrictions on the mean 
values of the streamwise velocity. 



For a time dependent solution, either the flux-vector-splitting or the flux- 
difference-splitting scheme is used in all three computational directions. However, 
these schemes are only used in T| and £ directions to obtain a space marching solution. 
In this case, either a PNS formulation or modified Roe's scheme is used to calculate the 
fluxes in the ^ direction. 

With alternate forward and backward relaxation sweeps, a time-dependent 
solution can be obtained for general flow problems. A one-sweep solution can be 
obtained for supersonic/subsonic mixing problems using either the PNS formulation 
or the modified Roe’s scheme. In these procedures, a local time-like iteration is used 
until the residual at the local plane is reduced four orders of magnitude. Space 
marching solutions require much less computational time than fully time-dependent 
solutions, and these can also be used as initial conditions to reduce the computational 
effort required by time-dependent solutions. 

5. Adaptive Grid Technique 

In the present report, the adaptive grid technique introduced by Pao and Abdol- 
Hamid^ is used to analyze the aerodynamic characteristics of shock containing single 
jets. This strategy is based upon the monitoring surface and equidistribution concept 
by Eiseman39‘41. Flow variables are used as the monitoring functions for grid 
adaptation in the computational domain of interest. Grid density is governed by the 
geometric properties of the monitoring functions. In this method, a template grid is 
first generated by a geometrical function which concentrates the grid points to a 
circular zone surrounding the jet plume while leaving an adequate number of grid 
points in the farfield computational domain. For the jet exit plane, this initial grid is 
adapted to axial velocity (representative of the shear layer) and the pressure 
(representative of the shock location). Adaptive grid for subsequent axial planes are 
generated by using the adapted grid in the previous plane as its template. 
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6. Boundary Conditions 


6.1 Multlblock/Multlzone Methodology 

The present multiblock/multizone strategy allows a great deal of topological 
flexibility. As long as there is no change in grid topology (polar or cartesian) or 
distribution (fine or coarse), additional blocks or zones are not necessary when the 
boundary condition changes along a block face. Each of the six block faces can have 
any combination of boundary conditions. Boundary conditions can be either a direct 
communication at the block interface or a regular boundary condition (inflow, 
outflow, reflection, symmetry wall or extrapolation). With this flexibility, the solver 
requires fewer blocks and zones which significantly reduces the overhead required for 
the communication between block/zone interfaces. 

A typical computational domain (figure 1) may contain zones, blocks and 
partitions. The relations between zones, blocks and partitions are defined as follows: 
Zone (Zj^): parts of the computational domain organized linearly along the i- 

direction. Each zone may contain multiple blocks. Each block within a zone can 
interface with blocks in adjacent zones where mixed boundary conditions are allowed. 
Communications at the zone level are restricted to face 5 (i zm = 1) which communicates 
with face 6 of zone Z m .i > and face 6 (i=i=imax zm ) which communicates with face 5 of 
the next zone 'Z m+ \. Blocks within a zone can terminate only at face 6 of the zone. Face 

5 of zone 1 and face 6 of the last zone of the computational domain can have any 
combination of the regular boundary conditions. 

Block (Bn): subsection of Zm, can start at any i-location (IMINzm < iBn < 
IMAXz m ) with different regular boundary conditions and interface with any other 
block (if needed) at faces 1.2.3 and 4. Face 5 and 6 can communicate with another from 
zone Z m _i and Z m+ i with any combination of regular boundary conditions. Face 1,2,3, 
and 4 are defined as follows: for i = Ign. imaxzm 
face 1 j=jmax, k=l, km ax 
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face 1 j=jmax, k=l, km ax 
face 2 j=jmin, k=l, km ax 
face 3 k=kmin, j=l, jmax 
face 4 k=kmax, j=l, jmax 

The interface between blocks can be of any direction or order (coarse to fine grid). Any 
of the four faces (1-4) in one block can communicate with any number of blocks and 
different faces in any order. With this flexibility, there is no need to break a block into 
a number of blocks on account of communications with more than one block or mixed 
boundary conditions on a block interface. The only computational effort is setting the 
correct boundary conditions for each of the four faces of a block. 

Partition (Pj): subsection of a block within a zone, which can start any i- 

location. The partition allows the user to change boundary conditions without adding 
an extra zone. Any of the four faces boundary conditions can be changed and a 
turbulent or laminar solution can be selected for each partition. 

Different examples of using the present multiblock/multizone methodology has been 
mentioned in Section 1. This method will simplify grid generation by reducing the 
number of blocks and zones needed to describe a complex computational domain and by 
giving the users more flexibility in breaking the computational domain into simple 
sections. 

6.2 Regular Boundary Conditions. 

Regular boundary conditions are inflow, outflow, reflection, symmetry, wall 
and extrapolation. In the following equations, a boundary point is denoted by "iy\ and 
the factitious image point is denoted by T. Boundary points (1) are calculated as 
follows: 

a) Inflow /Outflow Boundary 

The treatment of this boundary is based on Riemann invariants for a one- 
dimensional (low. Riemann Invariants can be defined as 
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( 41 ) 


f = v -i ± 7 L T°» 




where. + and - are the increasing and decreasing direction side 
of the ii-direction, and side of the h-direclion.. 

V,, = k x u + k y v + k z w 


then, there follows: 

v„=±i(p+6) 


a R 



(P-Q) 


The primitive flow variables can be written in a general fonn: 

{pf.Uf.Vf.Wf.pf} 1 ={c 1 .p 00 ,u„,v„, w^,c 2 p f } 1 +c 3 {0.k x ,k y ,k z .0} T 

where Cj, C2. C3 are given different values in various cases. 

1 ) For inflow boundary conditions where V R < 0 . 
c 3 =IV R -VJ 

If the constant entropy restriction is applied, then 


C] = 


( 2 \ 

a R 
2 
a_ 


1- 1 


_ a R 
*"2 = ~~~ 


otherwise, if the constant pressure restriction is applied, then 


Cj - 1 , c 2 — 


a^ 

Y 


2 ) For outflow boundary conditions where , Vr > o. 
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c 3 =IV R -V„l 


If the constant entropy restriction is applied, then 


c i - 


( 2 

a R 

2 

V a, i 


P3. 

Pr 


7-1 


,C 2 = 


a2 

R 

Y 


constant pressure restriction 

2 

a„, 

Y 

b) Reflection Boundary Condition 
There are three cases: 

u-reflection 

{p f , u f , v f , w f , p f } T = {p„. - u,,, v„. w n , pj T 
v-reflection 

{pf • u f . v f . w f , p f } T = {p„, u„, - v„, w,, p}‘ 

uv-reflection (quarter-plane polar grid) 

{Pf » Uf, v f , w f , Pf }" 1 ={p I ,.-u T 1 .-v n .w^.p tl } 1 

c) Symmetry Boundary Condition 

{pf . u f , v f , w f , p r } T = {p„, u n . v^,w rPl [ T 

d) Wall Boundary Condition 

At the wall, velocity normal to the wall surface is set to zero and no-slip conditions are 
imposed. 


- YP- - 
c i =- L 2 '.c 2 = 

a R 
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u =v = w = o 

The pressure on the wall surface is obtained by setting the pressure gradient to zero and 
adiabatic wall condition is employed for temperature. 

e) Extrapolation Boundary Condition 
There are three options: 
zero-order 

Qf= Qt) 

First-order 

Qf = 2 Qq + Qq . i 

Second-order 

Qf = 1.5Qq + 0.5 Qq . j 

7. Numerical Results and Discussions 

In this section, some of the computational options available in the PAB3-v2 
code are used to simulate different flow problems. Three test cases, which use the 
multiblock/multizone, adaptive grid, time-dependent, space marching, parabolized 
Navier-Stokes, and turbulence model capabilities are presented. The average 
computational time on Cray-2 computer was 70 p s /grid point. The first test case is a 

group of calculations for underexpanded and overexpanded supersonic jets issued into 
still air from round, square, and elliptic nozzles using polar and Cartesian grids. The 
second case is the simulation of high pressure supersonic flow issued into still air from 
a two 5-tabs nozzles. Finally, subsonic and supersonic flows past a nonaxisymmelric 
afterbody and nozzle with either a solid simulated plume or a supersonic jet exhaust are 
computed and compared with experimental data. 

7.1 Supersonic . let plume 

The present three-dimensional code is used to predict the shock-cell and flow 
characteristics of both underexpanded and overexpanded supersonic jets issued into 
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still air. The present predictions for round nozzles are qualitatively compared with the 
experimental data of Love et. al. 42 for underexpanded supersonic jets. Quantitative 
comparisons are made with the experimental data of Norum and Seiner 42 for Mach 2 
underexpanded and overexpanded supersonic jets issued from round nozzles. The 
present space marching scheme and adaptive grid are also used to simulate a supersonic 
underexpanded jet issued from both square and elliptic nozzles. 

Since details of the initial jet profile are not available, all cases are computed 
using a top-hat profile at the jet exit. Free-stream Mach number in the ambient air is 
assumed to have a value of .05. Space marching and time-dependent solutions 
presented in this section are calculated by using the van Leer flux-vector-splitting 
scheme. 

7.1.1 Polar Grid: Round Nozzle 

In order to evaluate the capabilities of the present code, some of the important 
characteristics of mildly underexpanded supersonic jets are computed and compared 
with PNS predictions as well as the experimental data of Love et. al. 42 . These are the 
characteristics of the first shock-cell: 

a) The location (i) of the intersection between the incident shock wave 
jet centerline or the intersection can be a Mach disk.. 

b) The location of the expansion wave reflection at 
jet boundary, W. 

Figure (2) shows a typical quarter plan polar grid used to calculate jet plume flow 
issued from a round nozzle. In figure (3), the computational results for i/D and W/D are 
presented as a function of pressure ratio P, / P a , and compared with the experimental 

data of Love et. al 42 - The magnitude of i and W increases with increasing pressure ratio 
for all three Mach numbers. Excellent agreement is achieved between the experimental 
data and computational results. In order to obtain PNS solutions, the highest 
applicable value of CFL number is approximately 5. By using the SMS technique, the 
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applicable values of CFL number can be as high as 30. Figure (4) shows a typical 
convergence history of the L-2 norm of the residual for SMS and PNS solutions. It is 
observed in this figure that the residual dropped by 10 orders of magnitude in less than 
100 iterations for the SMS solution. However, the PNS solution behaves similar to the 
SMS (in the first 10 iterations), then the convergence rate deteriorates substantially. In 
most cases, SMS solutions using PAB3D-v2 take less than 80 percent of the 
computational time required by PNS to achieve a similar convergence history. For the 
free stream region surrounding a Mach 1.5 or 2.0 jet, PNS solutions were unstable for a 
CFL value higher than to a value of 2 and a Mach number less than 0.3. It is not clear 
why a converged solution for overexpanded supersonic jet cases could not be obtained 
using PNS methodology. On the other hand, converged solutions were obtained for a 
wide range of Mach numbers and pressure ratio using the SMS. It was noted that the 
PNS procedure was very sensitive to the inlet condition at the jet exit and that a large 
discontinuity always caused a departure solution. 

These results indicate that SMS is a robust scheme which gives efficient and 
accurate solutions. These solutions can also be used as initial conditions for time- 
dependent Thin-Layer Navier-Stokes calculations which can reduce the number of 
iterations required for converged solutions^ 7 . SMS is to be used to calculate the test 
cases described in this section (section 7). 

Figure (5) shows the periodic structure of the Mach contours for a sonic jet with 
different pressure ratios up to 10 jet-radii downstream of the jet exit. Higher pressure 
ratios produced less shock cells within the same distance. The sonic flow at the exit is 
expanded to a supersonic flow with a higher Mach number. The number of shock-cells 
decreases from 6 to 3 and the fully expanded Mach number increases from 1. 1 to 2.0 as 
static pressure ratio is increased from 1.2 to 2.0. 


27 



7. 1.1.1 Underexpanded Mach 2.0 Jet: Pj/P a = 1.45 

The space marching prediction using the modified Roe's scheme in the 
streamwise direction is compared with a fully time-dependent solution in this section. 
The grid size for this case was 151x1 1x41 and the time-dependent solution converged 
after 500 global iterations, which is approximately 2100 CPU Secs. On the other hand, 
the space marching solution took less than 100 CPU Secs on the Cray 2. Figure (6) 
shows the centerline pressure variation predictions using both time -dependent and 
space marching methods. The space marching technique gives a solution comparable 
to the time-dependent solution in less than 5 percent of the computer time. All the first 
and second shock-cell characteristics (spacing and strength) are captured with the space 
marching scheme. 

Figure (7) shows a comparison of the predicted results using mixing length (ML) 
and mixing length with compressibility correction factor (ML-CC) turbulence models 
(SMS solution) with the measured streamwise pressure variation along the jet center 
line. The jet was operated at a pressure ratio of 1.45 corresponding to a fully expanded 
Mach number of 2.24 and was issued from convergent-divergent nozzle with a design 
Mach number of 2. 

The measured static pressure distributions indicate a decay in the shock 
structure strength which is due to the interaction of shocks with the growing mixing 
layer. The ML turbulence model significantly overpredicts the shock-cell decay. It was 
expected that the incompressible turbulence model would not give a good prediction for 
compressible flow problems. The predicted result improved with the use of the 
compressibility correction factor of Chuech et. al. 17 . The computational result agrees 
reasonably well with the experimental data up to 30 jet-radii downstream of the jet- 
exit. 
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7. 1.1.2 Overexpanded Mach 2 Jet: Pj/P^ = .75 

Figure (8) shows the present code prediction capability for an overexpanded flow 
case and for ML-CC and ML turbulent solutions, respectively. The jet was operated at a 
pressure ratio of 0.75 and was issued from a convergent-divergent nozzle with a design 
Mach number of 2. 

The code predicts the irregular centerline pressure data with remarkable 
precision. Both turbulence models give very similar pressure distributions and agree 
well with available data. This suggests that the mixing layer boundary does not reach 
the jet centerline. The code shows a sharp compression shock at the exit of jet as 
expected for this kind of flow. This increases the pressure downstream of the shock as 
shown in figure (8). 

7.1.2 Cartesian Grid 

A non-circular jet plume at olT-design operating conditions (over- or 
underexpanded) may contain a very complex internal shock cell structure and the shear 
layer cross section goes through a complex sequence of shape transformations in the 
developing region of the jet plume. In this section, solutions obtained by using both 
fixed and adapted grids for underexpanded supersonic jet issued from round, square, 
and elliptic nozzles are presented. Cartesian topology is used for the grids because it 
offers excellent grid mobility for the adaptive grid cases. Calculations are made using 
the space marching scheme in the stream wise direction and van Leer flux-vector- 
splitting scheme in the crossplane directions with laminar flow assumption. Initial 
shear layer thickness is assumed to be .05 and . 1 of the jet -radii for adaptive and fixed 
grid respectively. The test case is for an underexpanded supersonic jet with design Mach 
number of 2 and pressure ratio of 1.45. 

7. 1.2.1 Round Nozzle 

With only 36x36 grid points, the adaptive grid scheme has sufficient grid density 
to represent the circular nozzle shear layer as shown in figure (9). As the shear layer 
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and shock front change their location, the grid will follow these changes. Fixed and 
adaptive grids were used for single and multiblock solutions of an underexpanded 
supersonic jet (Pj/P a = 1.45). Figure (10) shows the centerline pressure distribution 

using a single-block solution with fixed and adaptive grid. The fixed grid solution is 
exactly the same as the results presented earlier using a polar grid topology. A 
secondary pressure rise at z/r e = 4 appears stronger in the adaptive grid calculation 

(agrees with experimental data), and there is a slight upstream shift of the shock-cell 
position from the fixed to the adapted grid results. In general, the adapted grid solution 
agrees better with the experimental data. Figure (11) shows how well the adapted grid 
follows the shear layer and detected shock front locations as clearly described by the 
density contours shown in figure (12). 

As can be noted in figure (9), cartesian grid topology wastes a large number of 
grid points in the farfleld region. To eliminate this problem, a multiblock strategy can 
be used in the farfield region as shown in figure (13). In this figure, the computational 
domain is divided into three blocks with the adaptive grid procedure applied only to the 
first block (with jet plume). Block dimensions are 27x27, 4x13, and 4x18 which results 
in about 33 percent less grid points and, in turn, causes a large reduction in 
computational time when compared to the single block strategy. As shown in figure 
(14), the three-block solution, which uses less computational resources (memory and 
time), agrees well with the result produced using the single block strategy. 

7. 1.2.2 Square and Elliptical Nozzle 

Figure (15) shows the initial adaptive grids generated for square and elliptic 
nozzle calculations. In the case of the square nozzle, the X-Z and Y-Z plane have 
similar shock-cell characteristics (density contours; figure (16)), to those presented for 
the round nozzle in figure (12). However, the cross section plane goes through a very 
complex transformation as the jet changes shape from square to round to diamond 
shaped Similar observations are made by Anderson and Barber^ for a supersonic 



rectangular nozzle. Figure (17) shows the comparison between the adaptive and fixed 
grid predictions of the centerline pressure. The shock-cell spacing and strength are 
very similar to the round nozzle results. Again, the secondary pressure strength at X/r e 

= 4 is much larger for the adaptive grid prediction than the fixed grid. In addition, there 
is another secondary pressure peak shown at X re = 1 7 for the adaptive grid. 

Figure (18) shows the density contours predictions for 2:1 elliptic nozzle. The 
shock-cell structure is very different from those presented for either the square or 
round nozzle. First, the minor axis grows each time that the shock reflected from the 
jet centerline intersects with the shear layer. At the same time, the major axis decays 
as the jet cross section transforms into a circular shape. Second, the shock front on the 
major axis plane (X-Z) is much stronger than the ones on the minor axis plane (Y-Z). 
Finally, the jet produces two different shock front structures associated with different 
scales (strength and space). It is believed that these structures have different centers 
which are not always located on the jet centerline (Z-axis). Figure (19) shows the 
centerline pressure distribution using fixed and adaptive grid solutions. The elliptic 
nozzle produces more shock-cells than either the round or square jet for a similar 
distance. In addition, the shock-cells have different centerline pressure characteristics 
(compare figure (19) with figures (10) and 17)). The adaptive grid solution predicts a 
large pressure peak for the first shock-cell whereas the fixed grid solution did not. 

7.2 Supersonic Nozzle 

In this section, the multiblock and partitioning options of PAB3D-v2 code are utilized 
to simulate a group of underexpanded supersonic jet flow which include both internal 
and external flow regions for a special family of jet nozzles. The examples are designed 
to show the flexibility of the PAB3D-v2 code in handling mixed boundary conditions 
over a block interface. The nozzle configuration can be described as a circular pipe 
section followed by five equally spaced tabs. Each tab is simply an extension of an arc 
segment of the circular pipe for a certain length in the downstream direction. Each arc 


segment, representing the width of the tab. is 1/10 of the full circle. Two nozzle 
configurations are chosen. In the first case, the tabs are infinite in length in a direction 
parallel to the pipe centerline, (fig. 20). This nozzle will be referred to as the "infinite 
tab nozzle". Essentially, the entire flow field can be considered as internal to the nozzle 
although there is venting through the spaces between the tabs. In the second case, the 
tabs are chosen to be one pipe diameter i length, (Fig. 21). The nozzle exit plane is 
defined as the streamwise location at the end of the tabs. This configuration is called 
the "short tab nozzle". The flow region upstream of the nozzle exit plane is considered 
as the nozzle internal flow, and the flow downstream of the same plane is a free jet with 
a shear layer developing between the jet plume and the ambient free stream. Each of the 
two cases has an internal (low Mach number of 2 with an internal to external pressure 
ratio of 1.45 at the nozzle exit. The cross sectional computational domain of the nozzle 
requires the extensive use of multiblock strategy. Most CFD codes with multiblock 
capability do not have the capability of handling mixed interface condition at a block 
boundary. Therefore, these codes would have to divide this computational domain into 
at least 15 blocks as shown in figure (22) with the number of blocks increasing with the 
number of tabs. However, PAB3D-v2 code needs to use only 2 blocks to simulate this 
test case for any number of tabs. For the short tab nozzle test case, most CFD codes need 
another block to simulate the flow downstream of the nozzle exit. Using the partition 
option, the PAB3D-v2 code does not have any need for an extra block or zone. With this 
flexibility, more than a 20 percent reduction in grid points is achieved: the overhead 
due to interpolation between block boundaries is avoided: and the convergence rate of 
the solution is increased. 

Figure (23) shows the Mach contours at different axial locations for the infinite 
tab nozzle test case. Initially, the flow escapes through the gaps between the 5 tabs and 
then the jet cross section changes into a pentagon shape connected at the midpoint of 
each tab. For, the short tab nozzle jet case (figure 24), the flow initially, goes through 
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the same transformation, but. the cross section changes its shape back to a circular one. 
Wlezien and Kibens 44 did a series of experiments for 1 ,2,4, and 8 tabs 
nonaxisymmetric nozzle-free jet. Their results indicate that multiple-tab nozzle 
plumes spread faster and have a lower core Mach number than the free-jet case for 
circular nozzles without tabs. The present predictions produce similar results to those 
reported in Ref. [44]. Figures (25) and (26) show the comparisons between centerline 
Mach number and pressure respectively for free-jet, nozzle-free jet, and nozzle test 
cases. The short tab nozzle case produces more shock-cells, decays faster, and has lower 
core Mach number than the free-jet described earlier in section 7 of this report. 

However, the infinite tab nozzle test case decays much faster, and has the lowest core 
Mach number of all cases considered. 

7.3 Nonaxisymmetric Afterbody 

In this section, the thin-layer Navier-Stokes equations are solved to simulate 
subsonic (M = 0.8) and supersonic (M = 1.2) flow cases over the nonaxisymmetric 
afterbody described in Ref. [45]. The third-order flux-difference splitting scheme of Roe 
is used to discretize the governing equations. In the subsonic case, jet-exhaust is 
included to utilize as an example to demonstrate the multiblock/multizone and 
partitioning options in the PAB3D-v2 code. The supersonic test case compares the 
performance of three different turbulence models in simulating a separated flow region. 

First, for an external flow Mach number of 0.8 and the jet operating at Mach 2, a 
coarse grid topology is used for the calculation (figure (27)). A laminar flow assumption 
is made for the regions upstream of the body and downstream of the jet exit. A turbulent 
flow calculation using the Baldwin-Lomax turbulence model (Ref. (10)) is made in the 
region over the body. Three different multiblock/multizone configurations (1Z-2B, 2Z- 
2B, 2Z-3B) are used to grid this problem. The first configuration has one zone with two 
blocks. The first block has 64x32x10 grid points to represent the full external flow 
region with three partitions ending at i = 10, 50, and 64. The second block simulates the 
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jet-exhaust region with 14x10x21 grid points. In the second configuration, two-zone 
block topology is used. The first zone has one block (50x10x32) with two partitions 
ending at i = 10 and 50. The second zone has one block (14x10x52) which represents the 
internal and external flow regions downstream of the jet exit. The last configuration is 
similar to the second one except that the second zone splits into two blocks; one for the 
external flow region (14x10x32) and the other for the jet-exit region (14x10x21). It 
should be mentioned that CFD codes which permit only one boundary condition per 
block interface would require a minimum of 4 blocks to simulate the present test case. 
With the PAB3D-v2 flexibility in dealing with mixed boundary conditions, the present 
case can be solved with as little as two blocks. 

Figure (28) shows the Mach contour predictions using the 1Z-2B topology (the 
other two topologies give exactly the same results). The smooth transition of the 
density contours between zone 1 and zone 2 is apparent in this figure. This transition is 
only possible because of the fully conservative nature of the zonal/block boundary 
scheme. The afterbody pressure distributions computed with the three different 
multiblock/multizone topologies are compared with experimental data in figure (29). 

All three calculations give exactly the same predictions (in less than 50 global 
iterations) and agree well with the experimental data. 

The supersonic Mach 1 .2 case is computed using three different turbulence 
models with a rectangular solid sling replacing the jet-exhaust. The three turbulence 
models are; the Baldwin-Lomax 1 ® model, the Johnson and King model *2-14 as 
extended to three-dimensional flows by Abid 1 ^ and Abid, et. al. and the Goldberg 
model 15-16 Calculations are made with a single-zone/single-block topology with fine 
grid distribution (129x66x33). 

The afterbody pressure distributions computed with the three different 
turbulence models (Baldwin-Lomax, Johnson-King, and Goldberg) are compared with 
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experimental data in figure (30). All three turbulence models predict the shock at 
approximately the correct axial location. The Johnson-King and Goldberg models, 
however, give shock locations slightly upstream which agree better with the data. 

Both the Johnson-King and Goldberg models give a much better prediction of the 
"pressure plateau" in the overcompression region. The Baldwin-Lomax model fails to 
predict the "pressure plateau". The Johnson-King model is presently the best in 
predicting the plateau. 

While the Johnson-King model may be slightly better in predicting the "pressure 
plateau," it had to be calibrated for this flow regime. For flows with massive 
separation, Johnson^ noticed that very large values of the ratio of the nonequilibrium 
maximum shear stresses were generated using the Johnson-King model. He suggested 
limiting this ratio should be limited to a maximum value of 3.0. For the present case, 
Abdol-Hamid and Compton** found that o < 2.5 produces one of the best calculations 
when compared with the experimental data. Abid et. al. used a limiting value of 4.0 
in predicting the pressure distributions for the ONERA M6 wing. For the present 
calculation, it was noted (not shown) that the convergence with a = 2.5 is slower than 
using the Baldwin-Lomax model. Again, Abid et. al. made similar observations. An 
advantage of both the Goldberg and Baldwin-Lomax models is that they do not need to 
be calibrated for this flow. 

8. Summary 

The PAB3D-v2 code and its application to a variety of aerodynamic test 
problems have been discussed in this technical report. The code solves the three- 
dimensional simplified Navier-Stokes equations using the strong conservation form of 
the finite volume formulations. It uses two different flux-splitting schemes: van Leer's 
flux-vector-splitting and Roe's flux-difference-splitting. Also, the modified Roe's 
scheme is used to efficiently solve the governing equations in the steady state mode 
with a single global sweep. Several of the PAB3D-v2 options (multiblock/multizone, 



adaptive grid, different turbulence models, and lime-dependent and space marching 
numerical schemes) were applied to three different flow examples. 

For the supersonic jet plume problem, the code uses its multiblock, adaptive 
grid, time-dependent and space marching strategies to predict the shock-cell structure 
for mildly underexpanded and overexpanded supersonic jet plumes. The jet is issued 
from round, square and 2: 1 elliptic nozzles. The multiblock option reduces the 
computational time by at least 20 percent . The space marching strategy predicts exactly 
the shock-cell structure and when compared to a fully time-dependent strategy, 
requires less than 5 percent of the computational time to get a converged solution. The 
adaptive grid option enhances the prediction of the flow solver as it compares better 
with the experimental data. 

For the vented supersonic nozzle lest cases, the PAB3D-v2 code provides a very 
efficient multiblock interface strategy to solve the 5-tabs nozzle. Many CFD codes 
require a great number of blocks per tab to solve these kinds of problems with mixed 
boundaries. Only two blocks are required by PAB3D-v2 to solve this problem, 
regardless of the number of tabs and length. 

The time-dependent option was utilized to simulate the flowfields around a 
nonaxisymmetric afterbody with external Mach numbers of 0.8 and 1.2. For the 0.8 
Mach number case, the jet-exhaust was simulated with different mulliblock/multizone 
topologies. The PAB3D-v2 code requires only two blocks to solve the jet-exhaust 
problem while many CFD codes need to break the flowfield into at least four blocks. 
Three different turbulence models were used and evaluated in solving the Mach 1.2 test 
case. The results of this study show that the Johnson-King and Goldberg turbulence 
models give a much better prediction of the shock location and pressure plateau in the 
separated region than the Baldwin-Lomax model. 

In general, the PAB3D-v2 code can be used to simulate flowfields for complex 
aerodynamic configurations. Obviously, a detailed validation study using the several 
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options in the code needs to be made. However, there is still a need for more options and 
code upgrading. Upgrading the turbulence models from the algebraic level to the two- 
equation, multi-scale or algebraic Reynolds Stress level is needed to deal with 3D 
mixing and general aircraft configurations. Time-dependent adaptive grid 
methodology is needed for simulating vortical aspects of 3D mixing. Lastly, finite-rate 
chemistry with multiple species and multiphase solvers needs to be added to the code 
for solving jet plume and combustion problems. 
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Figure 1. PAB3D-v2 Typical Computational Domain. 
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Experimental data; ref. 42 
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Figure 3. Comparison of Predicted (3D PNS and SMS) and Measured First Shock-Cell 
Characteristics for Supersonic Free Jets. 
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Figure 5. Predicted Mach Contours (3D SMS) for Undemnpanded Sonic Jet Operated at 
different Pressure Ratios. 
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Figure 6. Comparisons Between Time-dependent and Space Marching Solutions in 
Predicting the Centerline Pressure of Underexpanded Mach 2 Jets and 
Pj/P a = 1.45. 
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Experimental data 
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Figure 9. Initial Adaptive Grid for a Round Nozzle, Underexpanded Exit Pressure Ratio. 
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Figure 10. Comparison Belween Adaptive and Fixed Grid Calculations In Predicting 
the Centerline Pressure of Underexpanded Mach 2 Jet. 
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Figure 12. Three-Dimensional Density Contour Predictions for Underexpanded 
Mach 2 Round Nozzle Using Single Block Strategy. Pj/Pg = 1.45. 
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Figure 15. Initial Adaptive Grid for Square and Elliptic Nozzle. Underexpanded Exit 
Pressure Ratio. 
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Figure 17. Comparison Between Adaptive and Fixed Grid Calculations in Predicting 
the Centerline Pressure of Underexpanded Mach 2 Square Nozzle and 
P,/P a = 1.45 
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Figure 19. Comparison Delween Adaptive and Fixed Grid Calculations in Predicting 
the Centerline Pressure of Underexpanded Mach 2 Elliptic Nozzle and 
Pj/P a = 1.45. 
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Figure 20. Infinite Tab Nozzle Test Case Configuration. 
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Figure 24. Mach Contour Prediction for the Short Tab Nozzle at Different Axial Locations 
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Figure 25. Centerline Mach Number Comparisons Between Free-Jet, Infinite and Short Tab 
Nozzle Solutions. 
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Figure 26. Centerline Pressure Comparisons Between Frec-Jet. Infinite and Short Tab Nozzle 
oolutlons. 
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Figure 27. Three-Dimensional Computational Grid for Nonaxisymmetric Afterbody 
T est Case. 
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Figure 28. Density Contours Prediction for Mach 0.8 Case Using 1Z-2B Topology. 





Figure 29. Comparisons Between the Predictions of the Three Different 

Mulliblock/Moltlzone Topologies (1Z-2B, 2Z-2B and 2Z-3B) for Mach 0.8 
Test Case. 
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